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Abstract 

It  is  difficult  to  predict  stability  properties  of  a  finite  difference 
scheme.  It  has  to  be  investigated  through  the  roots  of  the  Z-trans- 
formed  and  Fourier  transformed  difference  scheme  (modal  equation). 
To  simultaneously  investigate  several  schemes  for  the  viscoelastic  wave 
equation,  it  is  possible  to  derive  the  modal  equation  with  parame¬ 
terized  coefficients.  Several  conditionally  stable  schemes  were  found, 
where  the  most  efficient  is  a  staggered  scheme  with  a  stability  condi¬ 
tion  closely  resembling  that  of  an  elastic  scheme. 


1  Introduction 

Finite  difference  schemes  are  often  used  to  numerically  solve  ordinary  and 
partial  differential  equations,  [7].  It  is  easy  to  create  such  a  scheme,  since 
derivatives  are  simply  exchanged  for  difference  ratios,  which  approximate  the 
derivatives  to  some  accuracy.  The  accuracy  is  found  through  Taylor  expan¬ 
sion.  There  is  no  standard  rule  for  how  these  schemes  are  created  and  thus 
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leaves  quite  some  room  for  experimentation.  The  accuracy  of  the  finite  differ¬ 
ence  scheme  is  easy  to  find  and  it  is  commonly  a  good  measure  for  how  well 
the  scheme  will  solve  the  differential  equation,  i.e. ,  the  higher  the  accuracy, 
the  closer  the  numerical  solution  will  be  to  the  true  solution.  For  time  depen¬ 
dent  problems  the  solution  at  earlier  time  is  used  to  update  the  solution  to 
the  next  time  level.  There  are  two  distinctive  techniques  to  update  the  solu¬ 
tion  in  time,  explicit  or  implicit  updating  (time  stepping).  For  implicit  time 
stepping  it  is  necessary  to  solve  an  equation  system  to  update  the  solution  in 
time.  Explicit  time  stepping  is  thus  naturally  faster.  Implicit  time  stepping 
tends  to  attenuate  higher  frequency  modes  in  the  solution  as  well.  Explicit 
time  stepping  is  always  unstable  (solution  grows  uncontrollably)  for  at  least 
some  choices  of  parameters,  e.g.,  time  step,  whereas  implicit  time  stepping 
is  stable  (for  well-posed  problems).  Explicit  time  stepping  is  preferred  for 
hyperbolic  problems,  since  it  does  not  necessarily  attenuate  the  solution  (hy¬ 
perbolic  systems  are  energy  preserving)  and  is  fast.  Hence,  it  is  necessary 
to  find  the  parameter  choices  for  which  the  explicit  scheme  is  not  unstable, 
i.e.  the  stability  condition.  This  can  prove  difficult  since  some  schemes  are 
unconditionally  unstable  (no  such  choices  exist)  and  only  some  conditionally 
stable.  The  problem  is  thus  to  find  a  conditionally  stable  scheme  and  the 
stability  condition. 

It  is  possible  to  investigate  the  stability  of  a  particular  finite  difference 
scheme  by  Fourier  transforming  the  scheme  in  space  and  Z-transforming  the 
scheme  in  time.  The  result  is  an  equation  which  depends  on  wavenumber. 
The  scheme  is  stable  if  the  absolute  value  of  the  roots  (poles  of  the  scheme)  of 
the  equation  are  all  equal  to  or  less  than  one  for  all  possible  wavenumbers  [7]. 
The  more  strict  formulation  states:  It  is  an  necessary  condition  for  stability 
or  It  is  a  sufficient  condition  for  instability  IF  the  absolute  value  of  one  or 
more  of  the  poles  is  greater  than  one.  The  method  implicitly  implies  that 
parameters  for  the  differential  equation  are  independent  of  space  (constants). 

Some  initial  work  to  solve  the  viscoelastic  wave  equation  with  an  explicit 
finite  difference  scheme  encountered  stability  problems,  [1],  [2].  To  find  a 
reliable  scheme  it  was  clearly  necessary  to  thoroughly  investigate  the  stability 
properties  of  several  different  realizations.  All  the  realizations  are  first  or 
second  order  accurate  in  time  and  fourth  order  accurate  in  space  and  based 
on  a  stress- velocity  formulation  with  one  memory  variable,  [1],  [9]. 

The  differential  equation  and  the  several  different  finite-difference  scheme 
realizations  will  be  covered  in  the  first  section  as  well  as  the  investigation 
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method.  The  following  section  presents  the  investigation  and  stability  results 
for  the  schemes.  Finally,  the  most  efficient  conditionally  stable  scheme  will 
be  thoroughly  analyzed  for  its  stability  condition. 

2  Viscoelastic  Wave  Propagation,  Discrete 
Schemes  and  Method  for  Investigation 

The  1-D  viscoelastic  wave  equation  for  a  constitutive  relation  corresponding 
to  a  standard  linear  solid  is, 

P,t  =  —K{  1  +  t)v<x  -  r 
r,t  =  — —  (r  +  I<  tv<x) 

Tv  ’ 

l 

V,t  = - P,x 

p 

[5].  It  is  an  initial  value  problem,  i.e.,  p  =  r  =  v  =  0,  t  <C  0.  K  is  the  bulk 
modulus,  p  is  the  density,  and  r  is  roughly  proportional  to  the  reciprocal  of  Q 
and  a  measure  of  attenuation,  [5].  ra  is  a  relaxation  time  and  determines  the 
frequency  range  where  there  is  most  attenuation,  p  is  the  pressure  (stress) 
and  v  is  the  particle  velocity,  r  is  the,  so  called,  memory  variable,  which 
reflects  the  ’memory’  of  the  viscoelastic  medium. 

Several  different  types  of  schemes  and  grids  can  be  used  to  discretize  the 
continuous  system,  Equation  (1).  The  general  scheme  created  here  can  be 
first  or  second  order  accurate  in  time  and  fourth  order  accurate  in  space. 
It  is  possible  to  fit  the  scheme  for  some  instances  of  parameter  values  on  a 
staggered  grid.  Acoustic/elastic  problems  have  generally  been  solved  with 
leap-frog  schemes  for  updating  in  time  [8].  Hence,  leap-frog  stencils  will  be 
used  for  to  update  p  and  v.  The  equation  for  the  memory  variable  r  might 
for  small  values  of  ra  be  stiff  and  it  is  natural  to  solve  the  equation  with  an 
Euler-Backward  or  Crank-Nicolson  scheme,  cf.,  [6].  With  these  specification 
a  fairly  general  finite  difference  scheme  for  Equation  (1)  can  be  constructed. 
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Pnm+1  =  Pnm~l-£K(  l+T)Do,xVl- 

2A t  (cic2r”+1  +  (1  -  ci)r"  +  Cl(l  -  c2)r”"1) 

(2)  C+1  =  C-/3-(H^(«+1  +  (l-a)C^)+  , 

^xDoAhb2v^(i  -  M<A(i  -  ^Kr1)) 

<+1  =  c1  -  -£-pDo,xpnm 

where  n  corresponds  to  time  level  (t  =  nAt)  and  m  to  space  (x  =  mAx). 
D0yX  is  a  discrete  centered  spatial  differentiation  operator.  The  parameters 
(3,  a,  b-i,  62,  cl5  c2  are  used  to  change  the  scheme,  f3  €  {  —  1,0},  a,  6l5  &2,  ci,  c2  € 
[0, 1],  The  parameters  do  not  affect  the  purely  acoustic  part  of  the  scheme, 
but  only  the  part  connected  to  the  memory  variable.  An  abundance  of  re¬ 
alizations  are  clearly  possible  even  though  the  possibilities  are  restricted 
through  the  small  parameter  space,  [3,  a  etc..  To  investigate  the  stability 
of  the  different  realizations,  the  scheme  [Equation  (2)]  is  Fourier  transformed 
in  space  and  ^-transformed  in  time,  creating  a  fifth  or  sixth  order  equation. 


z5+/3  qrip2T b\  b2C\  c2  -  z4q- f 

z*+0  [—2  -  ^2(1  +  r(l  +  <7r(M2(l  -  cx)  +  cxc2(l  -  b,))))]  + 

z2+f3  ^2Tqr^1  _  _  CiJ  +  ~  C2)  +  CxC2(l  -  62))]  + 

z2  [2 q  -  ?/>2(l  +  r)q]  +  z1+(3  [1  -  *l>2Tqr((  1  -  &x)ci(l  -  c2)  +  (1  -  ca)&x(l  -  62))]  + 
A  [-'i/’2T<3'A(l  -  62)cx(1  -  c2)\  -  q 
=  0, 

(3) 

where  ^  =  4V  ,  ,  =  Jrffl,  ,r  =  ,  =  (1  + 


/3)^,  \  —  Co^|.  c0  =  is  the  zero  frequency  velocity  and  k  is  the 

wavenumber.  if>  depends  on  the  spatial  accuracy  of  the  scheme  and  the  par¬ 


ticular  stencil.  If  the  absolute  value  of  the  solutions  to  Equation  (3)  is  less 
than  one,  then  the  corresponding  scheme  is  stable. 
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3  Stability  Investigation 

We  will  investigate  three  suites  of  different  schemes.  The  three  suites  are 
distinguished  by  the  values  of  f3  and  a.  These  parameters  changes  the  stencil 
for  updating  the  memory  variable  in  time.  The  three  different  combinations 
are:  (f3  =  0,a  =  1)  (Euler  backward),  (/?  =  0,  a  =  0.5)  (Crank-Nicolson), 
and  (/?  =  1  ,a  =  0.5)  (Crank-Nicolson). 

The  parameters  b{,Cj ;  i,j  6  {1,2},  determines  the  “mix”  in  time  of  the 
velocity  into  the  memory  variable  update  equation  and  of  the  memory  vari¬ 
able  into  the  pressure  update  equation. 

The  original  scheme  by  Blanch  et  ah,  [1],  [2],  was  equivalent  to  (3  =  0, 
a  =  0.5,  b\  =  0.5,  b2  =  1,  and  cx  =  0.  The  scheme  is  unconditionally 
unstable.  Figure  1  shows  the  position  of  the  poles  for  the  scheme  as  well 
as  other  related  schemes  with  /3  =  0,  which  results  in  a  fifth  order  equation 
and  five  poles.  A  first  attempt  to  adjust  the  scheme  might  be  to  change  the 
source  term  for  the  memory  variable  equation.  Changing  bx  between  0  and  1 
and  tracking  the  maximum  absolute  value  of  all  the  poles  yields  the  graph  in 
Figure  2.  Changing  a  to  a  =  1  to  create  an  even  more  stable  scheme  (Euler 
backward)  for  updating  the  memory  variable  in  time  has  no  effect  on  the 
stability  for  the  whole  scheme,  Figure  3.  Changing  cx  with  a  =  0.5,  6X  =  0.5, 
b2  =  1  ,  and  c2  =  0.5  produces  the  curve  in  Figure  4.  Figure  5  shows  the 
position  of  the  pole  which  causes  the  instability  of  the  scheme  for  small  values 
of  cx.  It  is  clearly  possible  to  create  a  conditionally  stable  scheme  for  the 
configuration  f3  =  0  and  a  —  0.5.  There  is  not  even  an  ad  hoc  explanation 
for  how  and/or  why  the  different  schemes  are  stable/unstable,  though.  Both 
bx  and  cx  can  be  varied  to  create  a  2-variable  plot  to  yield  more  insight  into 
the  phenomenom,  Figure  6.  Stable  schemes  occur  for  large  values  of  hx  and 
Ci,  suggesting  that  the  most  recent  information  should  be  used.  The  stable 
schemes  are  unfortunately  only  first  order  accurate,  though.  The  only  second 
order  accurate  scheme  occurs  for  bx  =  cx  =  0.5,  which  is  unstable. 

The  last  suite  of  schemes  (/?  =  1,  a  =  0.5)  are  for  some  instances  possible 
to  implement  as  staggered  schemes.  Fixing  b2  —  c2  =  0.5  will  ensure  that 
the  schemes  all  are  second  order  accurate  in  time  for  all  values  of  6X,  cx. 
The  positions  of  the  poles  for  this  suite  of  schemes  are  plotted  in  Figure 
7.  The  absolute  value  for  the  maximum  pole  as  a  function  of  6X  and  cx 
is  plotted  in  Figure  8  and  9.  Figure  10  shows  a  slice  from  Figure  8  for 
Ci  =  1.  Clearly  there  are  several  schemes  which  are  conditionally  stable. 
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Figure  1:  Position  of  poles  for  ft  =  0  schemes.  Unstable,  though  not  obvious. 
It  is  always  the  poles  close  to  —1,  which  are  outside  the  unit  circle  for  the 
unstable  schemes. 
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Figure  5:  The  position  of  the  maximum  pole  and  the  pole  corresponding 
to  the  largest  wavenumber.  They  coincide  for  small  values  of  c\.  For  large 
values  of  c\  the  largest  pole  corresponds  to  wavenumber  k  =  0  and  is  located 
at  1  +  0z\  Dots  and  crosses  conicides  until  the  scheme  is  stable.  Solid  -  Part 
of  the  unitcircle.  Crosses  -  The  position  of  the  largest  pole  (if  inside  scale  of 
plot).  Dots  -  Position  of  pole  corresponding  to  the  largest  wavenumber. 
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Figure  7:  Position  of  poles  for  /?  =  1  schemes,  for  a  large  value  of  the 
attenuation  parameter  r.  Clearly  unstable.  It  is  always  the  poles  close  to 
—  1,  which  are  outside  the  unit  circle  for  the  unstable  schemes. 

Two  schemes  (b\  =  0,ci  =  1)  and  (&i  =  l,ci  =  0)  are  possible  to  implement 
as  staggered  schemes.  One  with  the  memory  variable  r  located  at  the  same 
positions  as  the  pressure  p  or  located  at  the  same  positions  as  the  velocity 
v.  Figure  8  suggests  that  schemes  close  to  one  of  the  staggered  schemes  is 
conditionally  stable.  The  same  conclusion  could  be  made  from  the  pole  plot 
for  /?  =  0,  Figure  6,  even  though  it  is  not  as  clear  as  in  Figure  8,  since  the 
staggered  schemes  can  only  be  partially  constructed  for  /3  =  0. 

The  numerical  dispersion  properties  of  staggered  schemes  are  better  than 
non-staggered  schemes’,  which  make  a  staggered  formulation  more  prefferable 
than  a  non-staggered.  The  most  efficient  (least  necessary  amount  of  compu¬ 
tation)  staggered  formulation  places  the  memory  variable  r  at  the  location  of 
pressure  p.  This  is  the  most  efficient  scheme,  which  is  second  order  accurate 
in  time  and  fourth  order  accurate  in  space,  where  the  memory  variable  is 
updated  in  time  by  an  unconditionally  stable  scheme  (Crank-Nicolson). 
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Figure  8:  Absolute  value  of  maximum  pole  for  /?  =  1,  a  =  0.5,  62  =  C2  =  0.5 
and  varying  b{  and  C\.  Stable  schemes  occur  for  (equivalent)  staggered  sten¬ 
cils  and  schemes  with  coefficients  approaching  the  staggered  scheme  config¬ 
uration. 
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Figure  9:  Absolute  value  of  maximum  pole  for  (3  =  1,  a  =  0.5,  b2  =  c2  =  0.5 
and  varying  b\  and  C\.  Stable  schemes  occur  for  (equivalent)  staggered  sten¬ 
cils  and  schemes  with  coefficients  approaching  the  staggered  scheme  config¬ 
uration.  Other  view. 
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Figure  10:  Value  of  maximum  pole  minus  one.  Crank-Nicolson,  (3  =  1.  The 
scheme  becomes  unstable  for  b\  greater  than  approximately  0.47  <  b\. 
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4  Stability  Condition  for  the 
Staggered  Scheme 

It  is  necessary  to  know  the  stability  condition  of  a  finite  difference  scheme  in 
order  to  use  it  to  solve  a  differential  equation.  Hence,  the  stability  condition 
of  the  scheme  suggested  in  the  section  above  must  be  found.  The  Z-  and 
Fourier  transform  of  the  discrete  staggered  formulation  of  the  viscoelastic 
differential  equation  (1)  yields  the  following  third  order  equation,  [3],  [4], 
and  [9].  The  corresponding  sixth  order  equation,  Equation  (3),  contains  only 
even  order  terms  and  is  hence  equivalent  to  a  third  order  equation. 


z3  +  z2[-q  -  2  +  tp2(  1  +  r(l  -  qT))]  + 

(4)  z[l  +  2q-qip2(l  +  T(l  +  qr/q))\-q 

=  0 

where  q  =  ,r  =  fa,  end  4,  =  2A  ^  is  slightly 

different  from  before  due  to  the  staggered  grid.  There  are  a  few  simple 
limiting  cases  where  the  solution  easily  can  be  found.  If  rj  — >  0  then  Zit2  = 
1  -  ip2(  1  +  r)/2  ±  y^[l  -  ip2(l  +  r)/2]2  -  1,  z3  =  1,  and  for  k0  =  0,  z\y2  =  1 
and  Z3  =  q.  For  a  corresponding  acoustic  scheme  the  roots  would  be  zly2  = 

1  —  V,2/2  ±  —  V,2/2]2  —  1.  This  is  suggests  a  likely  stability  limit  at  either 

ip  <  2  or  ip\J  1  +  r  <  2.  The  limits  correspond  to 

(5)  a!  <  1 , 


where  A  could  be  either 

c0At 

Ax 
or 

/ f—j \  >  Cmax^t  Cq\/Y^VT At 

=  Ax  ~  Ax  ' 

The  number  cAt/ Ax  is  usually  referred  to  as  the  Courant  number. 

The  stability  limit  is  most  easily  investigated  numerically  by  fixing  the 
value  for  c0  and  Ax  and  varying  At.  Figure  11  shows  the  absolute  value  of 
the  maximum  pole  and  the  absolute  value  of  a  complex  pole  for  the  largest 
possible  wavenumber  (corresponding  to  spatial  Nyquist  frequency),  where 
t  —  0  and  T]  — >  0.  The  two  expressions,  Equation  (6)  and  (7),  for  A  yield 
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Figure  11:  Absolute  value  of  maximum  pole  and  absolute  value  of  pole  cor¬ 
responding  to  maximum  wavenumber.  The  poles  become  larger  than  one 
exactly  at  the  stability  limit,  r  =  0,  t]  — >  0.  Solid  -  Absolute  value  of  pole 
correspoding  to  maximum  wavenumber.  Plus  -  The  absolute  value  of  the 
pole  with  the  largest  absolute  value  of  all  poles  for  all  wavenumbers. 


the  same  stability  limit  for  r  =  0,  which  is  a  purely  acoustic  medium.  Using 
expression  (6)  for  the  stability  limit  for  r  =  0.02  ( Q  ~  100),  77  — +  0,  and 
plotting  the  same  poles  shows  that  the  acoustic  stability  limit  is  larger  than 
the  actual,  Figure  12.  Plotting  the  same  poles  but  using  the  expression  (7) 
as  the  stability  limit  results  in  Figure  13.  The  poles  become  larger  than  one 
exactly  at  the  stability  limit  in  Equation  (7). 

The  numerical  stability  results  supports  the  theory  for  77  — »•  0.  In  real 
simulations,  however,  77  will  be  a  finite  value,  since  otherwise  the  simulations 
will  essentially  be  acoustic  with  very  little  viscous  influences  unless  r  is  very 
large  (r  >  1).  Solving  Equation  (4)  for  77  —  0.2  and  r  =  0.02  yields  the  pole 
configuration  in  Figure  14.  The  stability  limit  seems  to  be  the  same  for  a 
finite  value  of  77.  Increasing  r  tor  =  0.2  has  no  effect  on  the  stability  limit,  yet 
is  the  configuration  changed,  Figure  15.  It  is  possible  to  see  in  Figure  15  that 
the  absolute  value  of  the  pole  corresponding  to  the  maximum  wavenumber 
starts  moving  towards  one  (the  unitcircle)  for  lower  wavenumbers  than  in 
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Courant  number  as  fraction  of  stability  criteria 

Figure  12:  Absolute  value  of  maximum  pole  and  absolute  value  of  pole  cor¬ 
responding  to  maximum  wavenumber.  The  Courant  number  is  compared 
to  the  stability  limit  in  Equation  (6).  The  scheme  becomes  unstable  for  a 
Courant  number  slightly  less  than  that  stability  limit,  r  =  0.02,  r)  —>  0. 
Solid  -  Absolute  value  of  pole  correspoding  to  maximum  wavenumber.  Plus 
-  The  absolute  value  of  the  pole  with  the  largest  absolute  value  of  all  poles 
for  all  wavenumbers. 
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Figure  13:  Absolute  value  of  maximum  pole  and  absolute  value  of  pole  cor¬ 
responding  to  maximum  wavenumber.  The  poles  become  larger  than  one 
exactly  at  the  stability  limit  in  Equation  (7).  r  =  0.02,  77  — >  0.  Solid  - 
Absolute  value  of  pole  correspoding  to  maximum  wavenumber.  Plus  -  The 
absolute  value  of  the  pole  with  the  largest  absolute  value  of  all  poles  for  all 
wavenumbers. 
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Figure  14:  Absolute  value  of  maximum  pole  and  absolute  value  of  pole  cor¬ 
responding  to  maximum  wavenumber.  The  poles  become  larger  than  one 
exactly  at  the  stability  limit  in  Equation  (7).  r  =  0.02,  rj  =  0.2.  Solid  - 
Absolute  value  of  pole  correspoding  to  maximum  wavenumber.  Plus  -  The 
absolute  value  of  the  pole  with  the  largest  absolute  value  of  all  poles  for  all 
wavenumbers. 

Figure  14.  For  r]  =  5  Figure  16  suggests  the  same  stability  limit  as  before. 
The  specific  behavior  of  the  poles  is  naturally  different  from  before. 

The  exact  pole  locations  for  a  range  of  wavenumbers  are  for  r  =  0.5, 
t]  =  0.2  shown  in  Figure  17  (stable)  and  18  (unstable).  The  complex  poles 
corresponding  to  propagating  modes  become  real  for  0  large  enough  to  violate 
the  stability  condition.  One  of  the  poles  attains  an  absolute  value  larger  than 
one  whereas  the  other  has  an  absolute  value  less  than  one.  Figures  19  and  20 
show  corresponding  plots  for  T]  =  5  (r  =  0.5).  Here  the  originally  real  pole 
breaks  out  of  the  unit  circle  in  contrast  to  the  i)  —  0.2  case. 
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Figure  15:  Absolute  value  of  maximum  pole  and  absolute  value  of  pole  cor¬ 
responding  to  maximum  wavenumber.  The  poles  become  larger  than  one 
exactly  at  the  stability  limit  in  Equation  (7).  r  =  0.2,  T)  =  0.2.  Solid  - 
Absolute  value  of  pole  correspoding  to  maximum  wavenumber.  Plus  -  The 
absolute  value  of  the  pole  with  the  largest  absolute  value  of  all  poles  for  all 
wavenumbers. 
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Figure  16:  Absolute  value  of  maximum  pole  and  absolute  value  of  pole  cor¬ 
responding  to  maximum  wavenumber.  The  poles  become  larger  than  one 
exactly  at  the  stability  limit  in  Equation  (7).  r  =  0.02,  77  =  5.  Solid  - 
Absolute  value  of  pole  correspoding  to  maximum  wavenumber.  Plus  -  The 
absolute  value  of  the  pole  with  the  largest  absolute  value  of  all  poles  for  all 
wavenumbers. 
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Figure  17:  Location  of  poles  (roots)  for  the  Courant  number  below  the  sta¬ 
bility  limit,  Tj  —  0.2  and  r  =  0.5.  Two  sets  of  poles  closely  follow  the  unit 
circle,  whereas  a  purely  real  pole  changes  value  around  q. 
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Figure  18:  Location  of  poles  (roots)  for  the  Courant  number  above  the  sta¬ 
bility  limit,  r]  =  0.2  and  r  =  0.5.  Two  sets  of  poles  closely  follow  the  unit 
circle,  whereas  a  purely  real  pole  changes  value  around  q.  The  poles  following 
the  unit  circle  pass  for  some  wavenumber  the  point  —  1  +  Oi  and  both  become 
purely  real.  One  with  an  absolute  value  greater  than  one  and  one  less  than 
one.  The  modes  corresponding  to  wavenumbers  where  the  poles  are  greater 
than  one  will  numerically  explode. 
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Real  part 

Figure  19:  Location  of  poles  (roots)  for  the  Courant  number  below  the  sta¬ 
bility  limit,  rj  =  5  and  r  =  0.5.  Two  sets  of  poles  closely  follow  the  unit 
circle  for  small  wavenumbers  but  have  considerably  smaller  absolute  value 
for  large  wavenumbers.  A  purely  real  pole  changes  value  around  q.  Here  q  is 
negative. 
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Real  part 


Figure  20:  Location  of  poles  (roots)  for  the  Courant  number  above  the  sta¬ 
bility  limit,  T]  =  5  and  r  =  0.5.  Two  sets  of  poles  closely  follow  the  unit 
circle  for  small  wavenumbers  but  have  considerably  smaller  absolute  value 
for  large  wavenumbers.  A  purely  real  pole  dcreases  its  value  from  q ,  which  is 
its  value  for  0  wavenumber.  The  purely  real  pole  passes  for  some  wavenum¬ 
ber  the  point  —1  +  0 i  and  attains  an  absolute  value  greater  than  one.  This  in 
contrast  to  the  case  t?  =  0.2  where  the  complex  poles  became  real  and  their 
absolute  greater  than  one.  The  modes  corresponding  to  wavenumbers  where 
the  pole  is  greater  than  one  will  numerically  explode. 
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5  Conclusions 


Stability  of  finite  difference  schemes  is  difficult  to  predict.  The  poles  of  the 
schemes  must  be  analyzed  for  all  possible  parameter  combinations.  Several 
schemes  can  be  investigated  simultaneously  by  deriving  the  modal  equation 
for  a  set  of  general  parameters. 

Several  conditionally  stable  finite  difference  scheme  were  found  for  the 
viscoelastic  wave  equation  corresponding  to  a  constitutive  relation  of  one 
standard  linear  solid.  The  most  numerically  efficient  scheme  is  a  staggered 
scheme  second  order  accurate  in  time  and  fourth  order  accurate  in  space. 
Only  spatially  fourth  order  accurate  schemes  were  investigated. 

The  stability  limit  for  the  staggered  scheme  is 

Cmax^i  ^  fi 
A-r  ~  1\[D  ’ 

where  D  is  the  dimension  and  cmax  is  the  infinite  frequency  velocity.  Different 
poles  break  out  of  the  unit  circle  for  too  large  Courant  numbers  depending 
on  the  value  77. 
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